Victorian Forest Futures Population Vability Analysis: Central Highlands

1 Project setup

1.0.1 Packages

Install packages if necessary. Listed are installers for packages not on CRAN and the devlopment (most up-to-date) version of STEPS.

# library(devtools)
# install_github("steps-dev/steps")
# install_github("smwindecker/gdaltools")

Load packages

library(raster)

Loading required package: sp

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:raster’:

    intersect, select, union

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(tibble)
library(sf)

Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3

library(rgdal)

rgdal: version: 1.4-3, (SVN revision 828)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
 Path to GDAL shared files: /usr/share/gdal/2.2
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.3-1 

library(readr)
library(readxl)
library(ggplot2)
library(lubridate)

Attaching package: ‘lubridate’

The following object is masked from ‘package:base’:

    date

library(magrittr)

Attaching package: ‘magrittr’

The following object is masked from ‘package:raster’:

    extract

library(tidyr)

Attaching package: ‘tidyr’

The following object is masked from ‘package:magrittr’:

    extract

The following object is masked from ‘package:raster’:

    extract

library(foreach)
library(doMC)

Loading required package: iterators
Loading required package: parallel

library(future)

Attaching package: ‘future’

The following object is masked from ‘package:raster’:

    values

library(future.apply)
library(tidyr)
library(dismo)
library(gbm)

Loaded gbm 2.1.5

library(steps)
library(gdaltools)

1.0.2 Functions

source(file = "R/functions/pg.sf.R")
source(file = "R/functions/pg.pa.R")
source(file = "R/functions/read.vba.R")
source(file = "R/functions/proc.vba.R")
source(file = "R/functions/get.landis.vars.R")
source(file = "R/functions/rascc.R")
source(file = "R/functions/read.multi.line.header.R")
source(file = "R/functions/interpolate.climdat.R")
source(file = "R/functions/get.rst.prop.R")
source(file = "R/functions/get.rst.dat.R")

1.0.3 Paths

proj_path <- "/home/landis/rfst/"
# proj_path <- "D:/Users/ryan/Dropbox/Work/RFA_STEPS/rfst/"

2 Prepare variables

2.1 Analysis parameters

2.1.1 Timesteps

ntimesteps <- 50

ncores <- 20

2.2 Geographic data

2.2.1 Landscape

Set up base landscape layer

ch_rst <- raster(x = "data/grids/eco_v12.img") %T>%
  plot

ch_proj <- ch_rst@crs
ch_extent <- extent(ch_rst)
ch_res <- res(ch_rst)

2.2.2 Boundaries

rfa <- read_sf("data/shapefiles/RFA/")%>%
  st_transform(crs = ch_proj)
rfa

Simple feature collection with 81 features and 7 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -57985.97 ymin: 5665387 xmax: 763088.7 ymax: 6223446
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

pg.sf(rfa)

ch_rfa <- rfa[rfa$NAME== "CENTRAL HIGHLANDS",] 
ch_rfa

Simple feature collection with 1 feature and 7 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 318074.3 ymin: 5770575 xmax: 452175.7 ymax: 5906852
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

pg.sf(ch_rfa)

ch_mask <- rasterize(ch_rfa, ch_rst, field = 1, filename = "./output/landscape_vars/ch_mask.grd", overwrite = TRUE) 
ch_mask[ch_rst == 132] <- NA # removes lake areas from habitat
plot(ch_mask)

Think carefully about using this layer as projections are for zone 55, so distorts west of the state

victoria <- read_sf("data/shapefiles/vicstatepolygon/") %>%
   st_transform(crs = ch_proj)

2.3 Occurrence data

pseudo_abs <- read_sf("data/shapefiles/pseudo_absences/pseudo_absences.shp")%>%
  st_transform(crs = st_crs(ch_rst)) %>%
  mutate(PA = 0) %>%
  select(PA, geometry) %>%
  st_zm(drop = TRUE)
pseudo_abs

Simple feature collection with 35 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 325312 ymin: 5774360 xmax: 445518.4 ymax: 5901431
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

2.3.1 GG Occurrence data

Read in data

gg_ari <- read_excel(path = "data/tabular/BoA_SB_Combined_VBA_upload_v4.xls") %>%
  dplyr::select(-starts_with("leave")) %>%
  rename("lon" = `X-coordinate (easting or longitude)`, "lat" = `Y-coordinate (northing or latitude)`) %>%
   fill(lon, lat, .direction = "down") %>%
  filter(`Taxon Name` == "Misc Target taxa not found") %>%
  select(lon, lat) %>%
  mutate(PA = 0) %>%
  st_as_sf(coords = c("lon", "lat"), crs = st_crs(28355)) %>%
  st_transform(crs = st_crs(ch_rst))

-
/
                                                                                                                  

New names:
* `leave blank` -> `leave blank...2`
* `leave blank` -> `leave blank...3`
* `leave blank` -> `leave blank...6`
* `leave blank` -> `leave blank...7`
* `leave blank` -> `leave blank...25`

gg_ari

Simple feature collection with 51 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 333214 ymin: 5779763 xmax: 451285 ymax: 5937215
epsg (SRID):    28355
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

gg_vba <- proc.vba("data/tabular/vba_gg_all_20190509.csv", project.crs = ch_proj)
gg_vba

Simple feature collection with 3086 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 248825.6 ymin: 5725915 xmax: 742929.5 ymax: 6008846
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

gg_pa_all <- gg_vba %>%
  rbind(gg_ari)
gg_pa_all

Simple feature collection with 3137 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 248825.6 ymin: 5725915 xmax: 742929.5 ymax: 6008846
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

plot_gg_pa_all <- pg.pa(gg_pa_all, ch_rfa) +
  geom_sf(data = victoria, fill = NA)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

plot_gg_pa_all

gg_pa_ch <- gg_pa_all[ch_rfa,]
gg_pa_ch

Simple feature collection with 1308 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 333053.7 ymin: 5783845 xmax: 448916.2 ymax: 5887474
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

gg_pa_ch_psab <- gg_pa_ch %>%
  rbind(pseudo_abs)

st_write(gg_pa_ch, "output/shapefiles/pa/ch/gg_pa_ch.shp", delete_layer = TRUE)

Deleting layer `gg_pa_ch' using driver `ESRI Shapefile'
Writing layer `gg_pa_ch' to data source `output/shapefiles/pa/ch/gg_pa_ch.shp' using driver `ESRI Shapefile'
features:       1308
fields:         1
geometry type:  Point

st_write(gg_pa_ch_psab, "output/shapefiles/pa/ch/gg_pa_ch_psab.shp", delete_layer = TRUE)

Deleting layer `gg_pa_ch_psab' using driver `ESRI Shapefile'
Writing layer `gg_pa_ch_psab' to data source `output/shapefiles/pa/ch/gg_pa_ch_psab.shp' using driver `ESRI Shapefile'
features:       1343
fields:         1
geometry type:  Point

pg_gg_pa_ch <- pg.pa(gg_pa_ch %>% mutate(pseudo_abs = "vba") %>% rbind(pseudo_abs %>% mutate(pseudo_abs = "pseudo_absences")), ch_rfa) +
  facet_grid(.~pseudo_abs)
pg_gg_pa_ch

2.3.2 LBP Occurrence data

lb_pa_all <- proc.vba("data/tabular/vba_lb_all_20190509.csv", project.crs = ch_proj)

plot_lb_pa_all <- pg.pa(lb_pa_all, ch_rfa) +
  geom_sf(data = victoria, fill = NA)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

plot_lb_pa_all

lb_pa_ch <- lb_pa_all[ch_rfa,]
lb_pa_ch

Simple feature collection with 1404 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 364618.3 ymin: 5798912 xmax: 447312.4 ymax: 5868055
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

lb_pa_ch_psab <- lb_pa_ch %>%
  rbind(pseudo_abs)

st_write(lb_pa_ch, "output/shapefiles/pa/ch/lb_pa_ch.shp", delete_layer = TRUE)

Deleting layer `lb_pa_ch' using driver `ESRI Shapefile'
Writing layer `lb_pa_ch' to data source `output/shapefiles/pa/ch/lb_pa_ch.shp' using driver `ESRI Shapefile'
features:       1404
fields:         1
geometry type:  Point

st_write(lb_pa_ch_psab, "output/shapefiles/pa/ch/lb_pa_ch_psab.shp", delete_layer = TRUE)

Deleting layer `lb_pa_ch_psab' using driver `ESRI Shapefile'
Writing layer `lb_pa_ch_psab' to data source `output/shapefiles/pa/ch/lb_pa_ch_psab.shp' using driver `ESRI Shapefile'
features:       1439
fields:         1
geometry type:  Point

pg_lb_pa_ch <- pg.pa(lb_pa_ch %>% mutate(pseudo_abs = "vba") %>% rbind(pseudo_abs %>% mutate(pseudo_abs = "pseudo_absences")), ch_rfa) +
  facet_grid(.~pseudo_abs)
pg_lb_pa_ch

save.image("output/session_images/landscape-occupancy.RData")

#load(file = "output/session_images/landscape-occupancy.RData")

2.4 Landis output

2.4.1 Scenario 1: Current planned burning, current timber harvesting, climate change under RCP 4.5, current bushfire regime (Buisness as usual)

s1_vars_landis <- get.landis.vars(
  scn_path = "~/landis_ch_s1_pb-th-cc_50/",
  proj_path = proj_path,
  out_path = "/output/habitat_vars/",
  scn_id = "s1",
  proj_mask = ch_mask,
  timesteps = ntimesteps,
  cores = ncores
  )

plot(s1_vars_landis[[1]])

plot(s1_vars_landis[[ntimesteps + 1]])

save(s1_vars_landis, file = "output/habitat_vars/s1_vars_landis")
# load(file = "output/habitat_vars/s1_vars_landis_fut")

2.4.2 Scenario 2: Current planned burning, current timber harvesting, climate change under RCP 8.5, increatred bushfire regime

2.4.3 Scenario 3: Current planned burning, current timber harvesting, no climate change, current bushfire regime

s3_vars_landis <- get.landis.vars(
  scn_path = "~/landis_ch_s3_pb-th-cc0_50/",
  proj_path = proj_path,
  out_path = "/output/habitat_vars/",
  scn_id = "s3",
  proj_mask = ch_mask,
  timesteps = ntimesteps,
  cores = ncores
  )

plot(s3_vars_landis[[1]])

plot(s3_vars_landis[[ntimesteps + 1]])

save(s3_vars_landis, file = "output/habitat_vars/s3_vars_landis")
# load(file = "output/habitat_vars/s3_vars_landis_fut")

2.4.4 Scenario 4: Current planned burning, no harvesting, climate change under RCP 4.5, current bushfire regime

s4_vars_landis <- get.landis.vars(
  scn_path = "~/landis_ch_s4_pb-th0-cc_50/",
  proj_path = proj_path,
  out_path = "/output/habitat_vars/",
  scn_id = "s4",
  proj_mask = ch_mask,
  timesteps = ntimesteps,
  cores = ncores,
  harvest_timber = FALSE
  )

plot(s4_vars_landis[[1]])

plot(s4_vars_landis[[ntimesteps + 1]])

save(s4_vars_landis, file = "output/habitat_vars/s4_vars_landis")
# load(file = "output/habitat_vars/s4_vars_landis_fut")

2.5 ARI Data

2.5.1 Anisotronic heating x ruggedness

ahr <-raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/Anisotrophic_Heating_Ruggedness") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_ahr.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.2 Log verticical distance all wetlands

lvdaw <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/log_vertical_distance_all_wetlands_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_lvdaw.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.3 Log vertical distance major streams

lvdma <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/log_vertical_distance_major_streams_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_lvdma.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.4 Log vertical distance minor streams

lvdmi <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/log_vertical_distance_minor_streams_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask) %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.5 Log vertical distance saline wetlands

lvdsw <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/log_vertical_distance_saline_wetlands_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_lvdsw.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.6 Relative annual insolation

rai <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/relative_annual_insolation_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_rai.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.7 Thorium

tho <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/sept2014thorium") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_tho.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.8 Visible sky index

vsky <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/visible_sky_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_vsky.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.9 Topographic wetness index

twi <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/wetness_index_topocrop_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_twi.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.10 Wind exposition

Problem with layer

# winex <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/wind_exposition2015") %>%
  # projectRaster(to = ch_mask) %>%
  # mask(mask = ch_mask, filename = "output/habitat_vars/ari_winex.grd") %T>%
  # plot

2.5.11 ARI variables stack

2.5.11.1 Initial variables

ari_iv <- stack(lvdaw, lvdma, lvdmi, lvdsw)
names(ari_iv) <- c("lvdaw", "lvdma", "lvdmi", "lvdsw")

2.5.11.2 Future variables

Repeates into list for each year

ari_v <- vector("list", ntimesteps + 1)
for(i in 1:(ntimesteps+1)){
  ari_v[[i]] <- ari_iv
}

2.6 Climactic data

pc ~ percentage change ac ~ absolute change

tmax ~ max temperature tmin ~ minimum temperature prec ~ precipitation evtr ~ evapotranspiration

01 ~ January 07 ~ July

4.5 ~ rcp 4.5 8.5 ~ rcp 8.5

2.6.1 Current climate

From WorldClim

2.6.1.1 Precipitation in January

prec01 <- raster("data/grids/wc2.0_30s_prec_01.tif") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/clim_vars/prec01.grd")

3 projected point(s) not finiteno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

names(prec01) <- "prec01"
plot(prec01)

2.6.1.2 Precipitation in July

prec07 <- raster("data/grids/wc2.0_30s_prec_07.tif") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/clim_vars/prec07.grd")

3 projected point(s) not finiteno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

names(prec07) <- "prec07"
plot(prec07)

2.6.1.3 Max temperature in January

tmax01 <- raster("data/grids/wc2.0_30s_tmax_01.tif") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/clim_vars/tmax01.grd")

3 projected point(s) not finiteno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

names(tmax01) <- "tmax01"
plot(tmax01)

2.6.1.4 Minimum temperature in July

tmin07 <- raster("data/grids/wc2.0_30s_tmin_07.tif") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/clim_vars/tmin07.grd")

3 projected point(s) not finiteno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

names(tmin07) <- "tmin07"
plot(tmin07)

2.6.2 Future climate

Data from Climate Change in Australia

2.6.2.1 Change data

CV has updated ##### Absolute change in temperature

#absolute change in temp
raw_tmax01_4.5_ac <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/tasmax_Amon_CSIRO-Mk3-6-0_rcp45_r1i1p1_abs-change-wrt-seasavg-clim_native.csv?tasmax_january[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

raw_tmax01_8.5_ac <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/tasmax_Amon_CSIRO-Mk3-6-0_rcp85_r1i1p1_abs-change-wrt-seasavg-clim_native.csv?tasmax_january[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

raw_tmin07_4.5_ac <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/tasmin_Amon_CSIRO-Mk3-6-0_rcp45_r1i1p1_abs-change-wrt-seasavg-clim_native.csv?tasmin_july[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]") 

raw_tmin07_8.5_ac <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/tasmin_Amon_CSIRO-Mk3-6-0_rcp85_r1i1p1_abs-change-wrt-seasavg-clim_native.csv?tasmin_july[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

Years data

n_tmax01_4.5 <- as.numeric(sub("-.*", "", unique(raw_tmax01_4.5_ac$time)))
n_tmax01_8.5 <- as.numeric(sub("-.*", "", unique(raw_tmax01_8.5_ac$time)))
n_tmin07_4.5 <- as.numeric(sub("-.*", "", unique(raw_tmin07_4.5_ac$time)))
n_tmin07_8.5 <- as.numeric(sub("-.*", "", unique(raw_tmin07_8.5_ac$time)))

tmax01_4.5_ac <- rascc(raw_tmax01_4.5_ac, new.proj.layer = ch_mask)
tmax01_8.5_ac <- rascc(raw_tmax01_8.5_ac, new.proj.layer = ch_mask)
tmin07_4.5_ac <- rascc(raw_tmin07_4.5_ac, new.proj.layer = ch_mask)
tmin07_8.5_ac <- rascc(raw_tmin07_8.5_ac, new.proj.layer = ch_mask)

2.6.3 Percentage change in precipitation

raw_prec01_4.5_pc <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/pr_Amon_CSIRO-Mk3-6-0_rcp45_r1i1p1_perc-change-wrt-seassum-clim_native.csv?pr_january[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

raw_prec01_8.5_pc <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/pr_Amon_CSIRO-Mk3-6-0_rcp85_r1i1p1_perc-change-wrt-seassum-clim_native.csv?pr_january[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

raw_prec07_4.5_pc <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/pr_Amon_CSIRO-Mk3-6-0_rcp45_r1i1p1_perc-change-wrt-seassum-clim_native.csv?pr_july[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

raw_prec07_8.5_pc <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/pr_Amon_CSIRO-Mk3-6-0_rcp85_r1i1p1_perc-change-wrt-seassum-clim_native.csv?pr_july[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

prec01_4.5_pc <- rascc(raw_prec01_4.5_pc, new.proj.layer = ch_mask)
prec01_8.5_pc <- rascc(raw_prec01_8.5_pc, new.proj.layer = ch_mask)
prec07_4.5_pc <- rascc(raw_prec07_4.5_pc, new.proj.layer = ch_mask)
prec07_8.5_pc <- rascc(raw_prec07_8.5_pc, new.proj.layer = ch_mask)

n_prec01_4.5 <- as.numeric(sub("-.*", "", unique(raw_prec01_4.5_pc$time)))
n_prec01_8.5 <- as.numeric(sub("-.*", "", unique(raw_prec01_4.5_pc$time)))
n_prec07_4.5 <- as.numeric(sub("-.*", "", unique(raw_prec07_4.5_pc$time)))
n_prec07_8.5 <- as.numeric(sub("-.*", "", unique(raw_prec07_4.5_pc$time)))

2.6.4 Absolute predicted values

2.6.4.1 Temperature

add new function that writes these climate layers to files.

2.6.4.1.1 Jan max temperature RCP 4.5

tmax01_4.5 <- tmax01_4.5_ac + tmax01
names(tmax01_4.5) <- names(tmax01_4.5_ac)
plot(tmax01_4.5)

2.6.4.1.2 Jan max temperature RCP 8.5

tmax01_8.5 <- tmax01_8.5_ac + tmax01
names(tmax01_8.5) <- names(tmax01_8.5_ac)
plot(tmax01_8.5)

2.6.4.1.3 July min temperature RCP 4.5

tmin07_4.5 <- tmin07_4.5_ac + tmax01
names(tmin07_4.5) <- names(tmin07_4.5_ac)
plot(tmin07_4.5)

2.6.4.1.4 July min temperature RCP 8.5

tmin07_8.5 <- tmin07_8.5_ac + tmax01
names(tmin07_8.5) <- names(tmin07_8.5_ac)
plot(tmin07_8.5)

2.6.4.2 Precipitation

2.6.4.3 January precipitation RCP 4.5

prec01_4.5 <- prec01 * (1 + prec01_4.5_pc/100)
names(prec01_4.5) <- names(prec01_4.5_pc)
plot(prec01_4.5)

2.6.4.4 January precipitation RCP 8.5

prec01_8.5 <- prec01 * (1 + prec01_8.5_pc/100)
names(prec01_8.5) <- names(prec01_8.5_pc)
plot(prec01_8.5)

2.6.4.5 July precipitation RCP 4.5

prec07_4.5 <- prec07 * (1 + prec07_4.5_pc/100)
names(prec07_4.5) <- names(prec07_4.5_pc)
plot(prec07_4.5)

2.6.4.6 July precipitation RCP 8.5

prec07_8.5 <- prec07 * (1 + prec07_8.5_pc/100)
names(prec07_8.5) <- names(prec07_8.5_pc)
plot(prec07_8.5)

2.6.5 Interploate prediction data

Currently just repeats predictions, need updating to weight predictions within 10 years, cf data from CSIRO

tmax01_4.5_int <- interpolate.climdat(tmax01, tmax01_4.5, ny = 50, nf = n_tmax01_4.5, n0 = 2019, varname = "tmax01")
tmax01_8.5_int <- interpolate.climdat(tmax01, tmax01_8.5, ny = 50, nf = n_tmax01_8.5, n0 = 2019, varname = "tmax01")
tmin07_4.5_int <- interpolate.climdat(tmin07, tmin07_4.5, ny = 50, nf = n_tmin07_4.5, n0 = 2019, varname = "tmin07")
tmin07_8.5_int <- interpolate.climdat(tmin07, tmin07_8.5, ny = 50, nf = n_tmin07_8.5, n0 = 2019, varname = "tmin07")

prec01_4.5_int <- interpolate.climdat(prec01, prec01_4.5, ny = 50, nf = n_prec01_4.5, n0 = 2019, varname = "prec01")
prec01_8.5_int <- interpolate.climdat(prec01, prec01_8.5, ny = 50, nf = n_prec01_8.5, n0 = 2019, varname = "prec01")
prec07_4.5_int <- interpolate.climdat(prec07, prec07_4.5, ny = 50, nf = n_prec07_4.5, n0 = 2019, varname = "prec07")
prec07_8.5_int <- interpolate.climdat(prec07, prec07_8.5, ny = 50, nf = n_prec07_8.5, n0 = 2019, varname = "prec07")

2.6.6 Climate variable sets

2.6.6.1 Initial climate variables

clim_iv <- stack(prec01, prec07, tmax01, tmin07)

2.6.6.2 Future climate variables

clim_fv_cc0 <- vector("list", 50)
for(i in 1:50){
  clim_fv_cc0[[i]] <- clim_iv
}


clim_fv_4.5 <- mapply(prec01_4.5_int, prec07_4.5_int, tmax01_4.5_int, tmin07_4.5_int, FUN = stack)
clim_fv_8.5 <- mapply(prec01_8.5_int, prec07_8.5_int, tmax01_8.5_int, tmin07_8.5_int, FUN = stack)

save.image(file = "output/input_variables.RData")

2.7 Distribution model variables

2.7.1 Species LANDIS variables

2.7.1.1 GG LANDIS vars

gglv <- c("ggf", "ggd", "hbt_1k")

2.7.1.2 LB LANDIS vars

lblv <- c("lbm", "hbt_3h")

2.7.2 Scenario 1

2.7.2.1 Greater Glider

2.7.2.1.1 Scenario 1 initial variables

s1_iv_gg <- stack(s1_vars_landis_init[[c(gglv)]], clim_iv)

save(s1_iv_gg, file = "output/habitat_vars/s1_iv_gg")

2.7.2.1.2 Scenaio 1 model data

s1_md_gg <- s1_iv_gg %>%
  raster::extract(y = st_coordinates(gg_pa_ch)) %>%
  cbind("PA" = gg_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s1_md_gg

save(s1_md_gg, file = "output/habitat_vars/s1_md_gg")

2.7.2.1.3 Scenario 1 future variables

s1_fv_gg <- s1_vars_landis_fut %>%
  sapply("[[", gglv) %>%
  mapply(clim_fv_4.5, FUN = stack)

save(s1_fv_gg, file = "output/habitat_vars/s1_fv_gg")

2.7.2.2 Leadbeater’s Possum

2.7.2.2.1 Scenario 1 initial variables

s1_iv_lb <- stack(s1_vars_landis_init[[c(lblv)]], clim_iv)

save(s1_iv_lb, file = "output/habitat_vars/s1_iv_lb")

2.7.2.2.2 Scenaio 1 model data

s1_md_lb <- s1_iv_lb %>%
  raster::extract(y = st_coordinates(lb_pa_ch)) %>%
  cbind("PA" = lb_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s1_md_lb

save(s1_md_lb, file = "output/habitat_vars/s1_md_lb")

2.7.2.2.3 Scenario 1 future variables

s1_fv_lb <- s1_vars_landis_fut %>%
  sapply("[[", lblv) %>%
  mapply(clim_fv_4.5, FUN = stack)

save(s1_fv_lb, file = "output/habitat_vars/s1_fv_lb")

2.7.3 Scenario 2

awaiting updated climate variables

2.7.4 Scenario 3

2.7.4.1 Greater Glider

2.7.4.1.1 Scenario 3 initial variables

s3_iv_gg <- stack(s3_vars_landis_init[[c(gglv)]], clim_iv)

save(s3_iv_gg, file = "output/habitat_vars/s3_iv_gg")

2.7.4.1.2 Scenaio 3 model data

s3_md_gg <- s3_iv_gg %>%
  raster::extract(y = st_coordinates(gg_pa_ch)) %>%
  cbind("PA" = gg_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s3_md_gg

save(s3_md_gg, file = "output/habitat_vars/s3_md_gg")

2.7.4.1.3 Scenario 3 future variables

s3_fv_gg <- s3_vars_landis_fut %>%
  sapply("[[", gglv) %>%
  mapply(clim_fv_cc0, FUN = stack)

save(s3_fv_gg, file = "output/habitat_vars/s3_fv_gg")

2.7.4.2 Leadbeater’s Possum

2.7.4.2.1 Scenario 3 initial variables

s3_iv_lb <- stack(s3_vars_landis_init[[c(lblv)]], clim_iv)

save(s3_iv_lb, file = "output/habitat_vars/s3_iv_lb")

2.7.4.2.2 Scenaio 3 model data

s3_md_lb <- s3_iv_lb %>%
  raster::extract(y = st_coordinates(lb_pa_ch)) %>%
  cbind("PA" = lb_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s3_md_lb

save(s3_md_lb, file = "output/habitat_vars/s3_md_lb")

2.7.4.2.3 Scenario 3 future variables

s3_fv_lb <- s3_vars_landis_fut %>%
  sapply("[[", lblv) %>%
  mapply(clim_fv_cc0, FUN = stack)

save(s3_fv_lb, file = "output/habitat_vars/s3_fv_lb")

2.7.5 Scenario 4

2.7.5.1 Greater Glider

2.7.5.1.1 Scenario 4 initial variables

s4_iv_gg <- stack(s4_vars_landis_init[[c(gglv)]], clim_iv)

save(s4_iv_gg, file = "output/habitat_vars/s4_iv_gg")

2.7.5.1.2 Scenaio 4 model data

s4_md_gg <- s4_iv_gg %>%
  raster::extract(y = st_coordinates(gg_pa_ch)) %>%
  cbind("PA" = gg_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s4_md_gg

save(s4_md_gg, file = "output/habitat_vars/s4_md_gg")

2.7.5.1.3 Scenario 4 future variables

s4_fv_gg <- s4_vars_landis_fut %>%
  sapply("[[", gglv) %>%
  mapply(clim_fv_4.5, FUN = stack)

save(s4_fv_gg, file = "output/habitat_vars/s4_fv_gg")

2.7.5.2 Leadbeater’s Possum

2.7.5.2.1 Scenario 4 initial variables

s4_iv_lb <- stack(s4_vars_landis_init[[c(lblv)]], clim_iv)

save(s4_iv_lb, file = "output/habitat_vars/s4_iv_lb")

2.7.5.2.2 Scenaio 4 model data

s4_md_lb <- s4_iv_lb %>%
  raster::extract(y = st_coordinates(lb_pa_ch)) %>%
  cbind("PA" = lb_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s4_md_lb

save(s4_md_lb, file = "output/habitat_vars/s4_md_lb")

2.7.5.2.3 Scenario 4 future variables

s4_fv_lb <- s4_vars_landis_fut %>%
  sapply("[[", lblv) %>%
  mapply(clim_fv_4.5, FUN = stack)

save(s4_fv_lb, file = "output/habitat_vars/s4_fv_lb")

3 Distribution model fit and prediction

3.1 Greater Glider

3.1.1 Scenario 1

s1_mod_gg <- gbm.step(data = s1_md_gg, gbm.x = 2:8, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s1_mod_gg)

#s1_ip_gg <- predict(s1_iv_gg, s1_mod_gg, type = "response", n.trees = s1_mod_gg$gbm.call$best.trees, filename = "output/habitat_pred/s1_fp_gg-000.tif")

#writeRaster(s1_ip_gg, "output/habitat_pred/s1_ip_gg.tif", overwrite = TRUE)
s1_ip_gg <- raster(x = "output/habitat_pred/s1_ip_gg.tif")

plot(s1_ip_gg)

registerDoMC(cores = 20)

s1_fp_gg <- foreach(i = seq_len(length(s1_fv_gg))) %dopar% {
  predict(s1_fv_gg[[i]],
          s1_mod_gg,
          type = "response",
          n.trees = s1_mod_gg$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s1_fp_gg-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s1_fp_gg <- stack(s1_fp_gg)

save(s1_fp_gg, file = "output/habitat_pred/s1_fp_gg")
writeRaster(s1_fp_gg, filename = "output/habitat_pred/s1_fp_gg.tif")

3.1.2 Scenario 3

s3_mod_gg <- gbm.step(data = s3_md_gg, gbm.x = 3:7, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s3_mod_gg)

 # s3_ip_gg <- predict(s3_iv_gg, s3_mod_gg, type = "response", n.trees = s3_mod_gg$gbm.call$best.trees, filename ="output/habitat_pred/s3_fp_gg-000.tif")

 # writeRaster(s3_ip_gg, "output/habitat_pred/s3_ip_gg.tif", overwrite = TRUE)

s3_ip_gg <- raster(x = "output/habitat_pred/s3_ip_gg.tif")

plot(s3_ip_gg)

registerDoMC(cores = 20)

s3_fp_gg <- foreach(i = seq_len(length(s3_fv_gg))) %dopar% {
  predict(s3_fv_gg[[i]],
          s3_mod_gg,
          type = "response",
          n.trees = s3_mod_gg$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s3_fp_gg-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s3_fp_gg <- stack(s3_fp_gg)

save(s3_fp_gg, file = "output/habitat_pred/s3_fp_gg")
writeRaster(s3_fp_gg, filename = "output/habitat_pred/s3_fp_gg.tif")

3.1.3 Scenario 4

s4_mod_gg <- gbm.step(data = s4_md_gg, gbm.x = 2:8, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s4_mod_gg)

# s4_ip_gg <- predict(s4_iv_gg, s4_mod_gg, type = "response", n.trees = s4_mod_gg$gbm.call$best.trees,  filename = "output/habitat_pred/4_fp_gg-000.tif")

# writeRaster(s4_ip_gg, "output/habitat_pred/s4_ip_gg.tif", overwrite = TRUE)

s4_ip_gg <- raster(x = "output/habitat_pred/s4_ip_gg.tif")

plot(s4_ip_gg)

registerDoMC(cores = 20)

s4_fp_gg <- foreach(i = seq_len(length(s4_fv_gg))) %dopar% {
  predict(s4_fv_gg[[i]],
          s4_mod_gg,
          type = "response",
          n.trees = s4_mod_gg$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s4_fp_gg-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s4_fp_gg <- stack(s4_fp_gg)

save(s4_fp_gg, file = "output/habitat_pred/s4_fp_gg")
writeRaster(s4_fp_gg, filename = "output/habitat_pred/s4_fp_gg.tif")

3.2 Leadbeater’s Possum

3.2.1 Scenario 1

s1_mod_lb <- gbm.step(data = s1_md_lb, gbm.x = 2:7, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s1_mod_lb)

# s1_ip_lb <- predict(s1_iv_lb, s1_mod_lb, type = "response", n.trees = s1_mod_lb$gbm.call$best.trees, filename = "output/habitat_pred/s1_fp_lb-000.tif")

# writeRaster(s1_ip_lb, "output/habitat_pred/s1_ip_lb.tif", overwrite = TRUE)

s1_ip_lb <- raster(x = "output/habitat_pred/s1_ip_lb.tif")

plot(s1_ip_lb)

registerDoMC(cores = 20)

s1_fp_lb <- foreach(i = seq_len(length(s1_fv_lb))) %dopar% {
  predict(s1_fv_lb[[i]],
          s1_mod_lb,
          type = "response",
          n.trees = s1_mod_lb$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s1_fp_lb-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s1_fp_lb <- stack(s1_fp_lb)

save(s1_fp_lb, file = "output/habitat_pred/s1_fp_lb")
writeRaster(s1_fp_lb, filename = "output/habitat_pred/s1_fp_lb.tif")

3.2.2 Scenario 3

s3_mod_lb <- gbm.step(data = s3_md_lb, gbm.x = 2:7, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s3_mod_lb)

# s3_ip_lb <- predict(s3_iv_lb, s3_mod_lb, type = "response", n.trees = s3_mod_lb$gbm.call$best.trees, filename = "output/habitat_pred/s3_fp_lb-000.tif")

# writeRaster(s3_ip_lb, "output/habitat_pred/s3_ip_lb.tif", overwrite = TRUE)

s3_ip_lb <- raster(x = "output/habitat_pred/s3_ip_lb.tif")

plot(s3_ip_lb)

registerDoMC(cores = 20)

s3_fp_lb <- foreach(i = seq_len(length(s3_fv_lb))) %dopar% {
  predict(s3_fv_lb[[i]],
          s3_mod_lb,
          type = "response",
          n.trees = s3_mod_lb$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s3_fp_lb-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s3_fp_lb <- stack(s3_fp_lb)

save(s3_fp_lb, file = "output/habitat_pred/s3_fp_lb")
writeRaster(s3_fp_lb, filename = "output/habitat_pred/s3_fp_lb.tif")

3.2.3 Scenario 4

s4_mod_lb <- gbm.step(data = s4_md_lb, gbm.x = 2:7, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s4_mod_lb)

# s4_ip_lb <- predict(s4_iv_lb, s4_mod_lb, type = "response", n.trees = s4_mod_lb$gbm.call$best.trees, filename = "output/habitat_pred/s4_fp_lb-000.tif")

# writeRaster(s4_ip_lb, "output/habitat_pred/s4_ip_lb.tif", overwrite = TRUE)

s4_ip_lb <- raster(x = "output/habitat_pred/s4_ip_lb.tif")

plot(s4_ip_lb)

registerDoMC(cores = 20)

s4_fp_lb <- foreach(i = seq_len(length(s4_fv_lb))) %dopar% {
  predict(s4_fv_lb[[i]],
          s4_mod_lb,
          type = "response",
          n.trees = s4_mod_lb$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s4_fp_lb-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s4_fp_lb <- stack(s4_fp_lb)

save(s4_fp_lb, file = "output/habitat_pred/s4_fp_lb")
writeRaster(s4_fp_lb, filename = "output/habitat_pred/s4_fp_lb.tif")

4 Population viability analysis

4.1 Greater Glider

gg_trans_mat <- matrix(c(0.00,0.00,0.25,
                         0.50,0.00,0.00,
                         0.00,0.85,0.85),
                       nrow = 3, ncol = 3, byrow = TRUE) # based on Possingham et al 1994
colnames(gg_trans_mat) <- rownames(gg_trans_mat) <- c('Newborn','Juvenile','Adult')

gg_stable_states <- abs( eigen(gg_trans_mat)$vectors[,1] / base::sum(eigen(gg_trans_mat)$vectors[,1]) ) 

4.1.1 Scenario 1

s1_hs_gg <- stack(s1_ip_gg, s1_fp_gg)

for(i in 1:51){
  s1_hs_gg[[i]][is.na(s1_hs_gg[[i]][])] <- 0
}

s1_hs_gg <- mask(s1_hs_gg, mask = ch_mask)

s1_hab_k_gg <- calc(s1_hs_gg[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s1_hab_k_gg) <- "carryingCapacity"

s1_gg_popN <- stack(replicate(ncol(gg_trans_mat), s1_hab_k_gg))

s1_gg_popN <- s1_gg_popN*gg_stable_states

s1_gg_idx <- which(!is.na(s1_hs_gg[[1]][]) & s1_hs_gg[[1]][] < 0.95)

s1_gg_pop <- s1_gg_popN
s1_gg_pop[!is.na(s1_gg_pop)] <- 0
s1_gg_pop[[1]][sample(s1_gg_idx, 10000)] <- 1
s1_gg_pop[[2]][sample(s1_gg_idx, 10000)] <- 1
s1_gg_pop[[3]][sample(s1_gg_idx, 10000)] <- 1

s1_gg_pop <- s1_gg_pop*ch_mask

names(s1_gg_pop) <- colnames(gg_trans_mat)

s1_gg_TotpopN <- sum(cellStats(s1_gg_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s1_gg_init_pop_size <- sum(s1_gg_pop)

s1_gg_landscape <- landscape(population = s1_gg_pop,
                          suitability = s1_hs_gg,
                          carrying_capacity = s1_hab_k_gg)

s1_gg_pop_dynamics <- population_dynamics(change = growth(transition_matrix = gg_trans_mat,
                                                       global_stochasticity = 0.1),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 8000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s1_gg_results_1_50 <- simulation(landscape = s1_gg_landscape,
                         population_dynamics = s1_gg_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s1_gg_results_1_50, file = "output/pva_results/s1_gg_results_1_50.Rds")

plot(s1_gg_results_1_50)

plot(s1_gg_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))

4.1.2 Scenario 3

s3_hs_gg <- stack(s3_ip_gg, s3_fp_gg)

for(i in 1:51){
  s3_hs_gg[[i]][is.na(s3_hs_gg[[i]][])] <- 0
}

s3_hs_gg <- mask(s3_hs_gg, mask = ch_mask)

s3_hab_k_gg <- calc(s3_hs_gg[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s3_hab_k_gg) <- "carryingCapacity"

s3_gg_popN <- stack(replicate(ncol(gg_trans_mat), s3_hab_k_gg))

s3_gg_popN <- s3_gg_popN*gg_stable_states

s3_gg_idx <- which(!is.na(s3_hs_gg[[1]][]) & s3_hs_gg[[1]][] < 0.95)

s3_gg_pop <- s3_gg_popN
s3_gg_pop[!is.na(s3_gg_pop)] <- 0
s3_gg_pop[[1]][sample(s3_gg_idx, 10000)] <- 1
s3_gg_pop[[2]][sample(s3_gg_idx, 10000)] <- 1
s3_gg_pop[[3]][sample(s3_gg_idx, 10000)] <- 1

s3_gg_pop <- s3_gg_pop*ch_mask

names(s3_gg_pop) <- colnames(gg_trans_mat)

s3_gg_TotpopN <- sum(cellStats(s3_gg_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s3_gg_init_pop_size <- sum(s3_gg_pop)

s3_gg_landscape <- landscape(population = s3_gg_pop,
                          suitability = s3_hs_gg,
                          carrying_capacity = s3_hab_k_gg)

s3_gg_pop_dynamics <- population_dynamics(change = growth(transition_matrix = gg_trans_mat,
                                                       global_stochasticity = 0.1),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 8000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s3_gg_results_1_50 <- simulation(landscape = s3_gg_landscape,
                         population_dynamics = s3_gg_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s3_gg_results_1_50, file = "output/pva_results/s3_gg_results_1_50.Rds")

plot(s3_gg_results_1_50)

plot(s3_gg_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))

4.1.3 Scenario 4

s4_hs_gg <- stack(s4_ip_gg, s4_fp_gg)

for(i in 1:51){
  s4_hs_gg[[i]][is.na(s4_hs_gg[[i]][])] <- 0
}

s4_hs_gg <- mask(s4_hs_gg, mask = ch_mask)

s4_hab_k_gg <- calc(s4_hs_gg[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s4_hab_k_gg) <- "carryingCapacity"

s4_gg_popN <- stack(replicate(ncol(gg_trans_mat), s4_hab_k_gg))

s4_gg_popN <- s4_gg_popN*gg_stable_states

s4_gg_idx <- which(!is.na(s4_hs_gg[[1]][]) & s4_hs_gg[[1]][] < 0.95)

s4_gg_pop <- s4_gg_popN
s4_gg_pop[!is.na(s4_gg_pop)] <- 0
s4_gg_pop[[1]][sample(s4_gg_idx, 10000)] <- 1
s4_gg_pop[[2]][sample(s4_gg_idx, 10000)] <- 1
s4_gg_pop[[3]][sample(s4_gg_idx, 10000)] <- 1

s4_gg_pop <- s4_gg_pop*ch_mask

names(s4_gg_pop) <- colnames(gg_trans_mat)

s4_gg_TotpopN <- sum(cellStats(s4_gg_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s4_gg_init_pop_size <- sum(s4_gg_pop)

s4_gg_landscape <- landscape(population = s4_gg_pop,
                          suitability = s4_hs_gg,
                          carrying_capacity = s4_hab_k_gg)

s4_gg_pop_dynamics <- population_dynamics(change = growth(transition_matrix = gg_trans_mat,
                                                       global_stochasticity = 0.1),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 8000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s4_gg_results_1_50 <- simulation(landscape = s4_gg_landscape,
                         population_dynamics = s4_gg_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s4_gg_results_1_50, file = "output/pva_results/s4_gg_results_1_50.Rds")

plot(s4_gg_results_1_50)

plot(s4_gg_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))

4.2 Leadbeater’s Possum

lb_trans_mat <- matrix(c(0.00, 0.00, 1.10,
                         0.59, 0.00, 0.00,
                         0.00, 0.59, 0.79),
                       nrow = 3, ncol = 3, byrow = TRUE) 
colnames(lb_trans_mat) <- rownames(lb_trans_mat) <- c('Newborn','Juvenile','Adult')

lb_stable_states <- abs( eigen(lb_trans_mat)$vectors[,1] / base::sum(eigen(lb_trans_mat)$vectors[,1]) ) 

4.2.1 Scenario 1

s1_hs_lb <- stack(s1_ip_lb, s1_fp_lb)

for(i in 1:51){
  s1_hs_lb[[i]][is.na(s1_hs_lb[[i]][])] <- 0
}

s1_hs_lb <- mask(s1_hs_lb, mask = ch_mask)

s1_hab_k_lb <- calc(s1_hs_lb[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s1_hab_k_lb) <- "carryingCapacity"

s1_lb_popN <- stack(replicate(ncol(lb_trans_mat), s1_hab_k_lb))

s1_lb_popN <- s1_lb_popN*lb_stable_states

s1_lb_idx <- which(!is.na(s1_hs_lb[[1]][]) & s1_hs_lb[[1]][] < 0.95)

s1_lb_pop <- s1_lb_popN
s1_lb_pop[!is.na(s1_lb_pop)] <- 0
s1_lb_pop[[1]][sample(s1_lb_idx, 3000)] <- 1
s1_lb_pop[[2]][sample(s1_lb_idx, 3000)] <- 1
s1_lb_pop[[3]][sample(s1_lb_idx, 3000)] <- 1

s1_lb_pop <- s1_lb_pop*ch_mask

names(s1_lb_pop) <- colnames(lb_trans_mat)

s1_lb_TotpopN <- sum(cellStats(s1_lb_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s1_lb_init_pop_size <- sum(s1_lb_pop)

s1_lb_landscape <- landscape(population = s1_lb_pop,
                          suitability = s1_hs_lb,
                          carrying_capacity = s1_hab_k_lb)

s1_lb_pop_dynamics <- population_dynamics(change = growth(transition_matrix = lb_trans_mat,
                                                       global_stochasticity = 0.3),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 2000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s1_lb_results_1_50 <- simulation(landscape = s1_lb_landscape,
                         population_dynamics = s1_lb_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s1_lb_results_1_50, file = "output/pva_results/s1_lb_results_1_50.Rds")

plot(s1_lb_results_1_50)

plot(s1_lb_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))

4.2.2 Scenario 3

s3_hs_lb <- stack(s3_ip_lb, s3_fp_lb)

for(i in 1:51){
  s3_hs_lb[[i]][is.na(s3_hs_lb[[i]][])] <- 0
}

s3_hs_lb <- mask(s3_hs_lb, mask = ch_mask)

s3_hab_k_lb <- calc(s3_hs_lb[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s3_hab_k_lb) <- "carryingCapacity"

s3_lb_popN <- stack(replicate(ncol(lb_trans_mat), s3_hab_k_lb))

s3_lb_popN <- s3_lb_popN*lb_stable_states

s3_lb_idx <- which(!is.na(s3_hs_lb[[1]][]) & s3_hs_lb[[1]][] < 0.95)

s3_lb_pop <- s3_lb_popN
s3_lb_pop[!is.na(s3_lb_pop)] <- 0
s3_lb_pop[[1]][sample(s3_lb_idx, 3000)] <- 1
s3_lb_pop[[2]][sample(s3_lb_idx, 3000)] <- 1
s3_lb_pop[[3]][sample(s3_lb_idx, 3000)] <- 1

s3_lb_pop <- s3_lb_pop*ch_mask

names(s3_lb_pop) <- colnames(lb_trans_mat)

s3_lb_TotpopN <- sum(cellStats(s3_lb_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s3_lb_init_pop_size <- sum(s3_lb_pop)

s3_lb_landscape <- landscape(population = s3_lb_pop,
                          suitability = s3_hs_lb,
                          carrying_capacity = s3_hab_k_lb)

s3_lb_pop_dynamics <- population_dynamics(change = growth(transition_matrix = lb_trans_mat,
                                                       global_stochasticity = 0.1),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 8000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s3_lb_results_1_50 <- simulation(landscape = s3_lb_landscape,
                         population_dynamics = s3_lb_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s3_lb_results_1_50, file = "output/pva_results/s3_lb_results_1_50.Rds")

plot(s3_lb_results_1_50)

plot(s3_lb_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))

4.2.3 Scenario 4

s4_hs_lb <- stack(s4_ip_lb, s4_fp_lb)

for(i in 1:51){
  s4_hs_lb[[i]][is.na(s4_hs_lb[[i]][])] <- 0
}

s4_hs_lb <- mask(s4_hs_lb, mask = ch_mask)

s4_hab_k_lb <- calc(s4_hs_lb[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s4_hab_k_lb) <- "carryingCapacity"

s4_lb_popN <- stack(replicate(ncol(lb_trans_mat), s4_hab_k_lb))

s4_lb_popN <- s4_lb_popN*lb_stable_states

s4_lb_idx <- which(!is.na(s4_hs_lb[[1]][]) & s4_hs_lb[[1]][] < 0.95)

s4_lb_pop <- s4_lb_popN
s4_lb_pop[!is.na(s4_lb_pop)] <- 0
s4_lb_pop[[1]][sample(s4_lb_idx, 3000)] <- 1
s4_lb_pop[[2]][sample(s4_lb_idx, 3000)] <- 1
s4_lb_pop[[3]][sample(s4_lb_idx, 3000)] <- 1

s4_lb_pop <- s4_lb_pop*ch_mask

names(s4_lb_pop) <- colnames(lb_trans_mat)

s4_lb_TotpopN <- sum(cellStats(s4_lb_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s4_lb_init_pop_size <- sum(s4_lb_pop)

s4_lb_landscape <- landscape(population = s4_lb_pop,
                          suitability = s4_hs_lb,
                          carrying_capacity = s4_hab_k_lb)

s4_lb_pop_dynamics <- population_dynamics(change = growth(transition_matrix = lb_trans_mat,
                                                       global_stochasticity = 0.1),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 8000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s4_lb_results_1_50 <- simulation(landscape = s4_lb_landscape,
                         population_dynamics = s4_lb_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s4_lb_results_1_50, file = "output/pva_results/s4_lb_results_1_50.Rds")

plot(s4_lb_results_1_50)

plot(s4_lb_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))